home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-02 / nrpas13.zip / AMOEBA.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  3KB  |  91 lines
  1. PROCEDURE amoeba(VAR p: glmpnp; VAR y: glmp; ndim: integer;
  2.        ftol: real; VAR iter: integer);
  3. (* Programs using routine AMOEBA must supply an external function
  4. func(pr:glnp):real whose minimum is to be found. They must
  5. also define types
  6. TYPE
  7.    glmpnp = ARRAY [1..mp,1..np] OF real;
  8.    glmp = ARRAY [1..mp] OF real;
  9.    glnp = ARRAY [1..np] OF real;
  10. where mp and np are physical dimensions *)
  11. LABEL 99;
  12. CONST
  13.    alpha=1.0;
  14.    beta=0.5;
  15.    gamma=2.0;
  16.    itmax=500;
  17. VAR
  18.    mpts,j,inhi,ilo,ihi,i: integer;
  19.    yprr,ypr,rtol: real;
  20.    pr,prr,pbar: glnp;
  21. BEGIN
  22.    mpts := ndim+1;
  23.    iter := 0;
  24.    WHILE true DO BEGIN
  25.       ilo := 1;
  26.       IF (y[1] > y[2]) THEN BEGIN
  27.          ihi := 1;
  28.          inhi := 2
  29.       END ELSE BEGIN
  30.          ihi := 2;
  31.          inhi := 1
  32.       END;
  33.       FOR i := 1 TO mpts DO BEGIN
  34.          IF (y[i] < y[ilo]) THEN  ilo := i;
  35.          IF (y[i] > y[ihi]) THEN BEGIN
  36.             inhi := ihi;
  37.             ihi := i
  38.          END ELSE IF  (y[i] > y[inhi]) THEN
  39.             IF (i <> ihi) THEN  inhi := i
  40.       END;
  41.       rtol := 2.0*abs(y[ihi]-y[ilo])/(abs(y[ihi])+abs(y[ilo]));
  42.       IF (rtol < ftol) THEN GOTO 99;
  43.       IF (iter = itmax) THEN BEGIN
  44.          writeln('pause in AMOEBA - too many iterations'); readln
  45.       END;
  46.       iter := iter+1;
  47.       FOR j := 1 TO ndim DO pbar[j] := 0.0;
  48.       FOR i := 1 TO mpts DO
  49.          IF (i <> ihi) THEN FOR j := 1 TO ndim DO pbar[j] := pbar[j]+p[i,j];
  50.       FOR j := 1 TO ndim DO BEGIN
  51.          pbar[j] := pbar[j]/ndim;
  52.          pr[j] := (1.0+alpha)*pbar[j]-alpha*p[ihi,j]
  53.       END;
  54.       ypr := func(pr);
  55.       IF (ypr <= y[ilo]) THEN BEGIN
  56.          FOR j := 1 TO ndim DO prr[j] := gamma*pr[j]+(1.0-gamma)*pbar[j];
  57.          yprr := func(prr);
  58.          IF (yprr < y[ilo]) THEN BEGIN
  59.             FOR j := 1 TO ndim DO p[ihi,j] := prr[j];
  60.             y[ihi] := yprr
  61.          END ELSE BEGIN
  62.             FOR j := 1 TO ndim DO p[ihi,j] := pr[j];
  63.             y[ihi] := ypr
  64.          END
  65.       END ELSE IF  (ypr >= y[inhi]) THEN BEGIN
  66.          IF (ypr < y[ihi]) THEN BEGIN
  67.             FOR j := 1 TO ndim DO p[ihi,j] := pr[j];
  68.             y[ihi] := ypr
  69.          END;
  70.          FOR j := 1 TO ndim DO prr[j] := beta*p[ihi,j]+(1.0-beta)*pbar[j];
  71.          yprr := func(prr);
  72.          IF (yprr < y[ihi]) THEN BEGIN
  73.             FOR j := 1 TO ndim DO p[ihi,j] := prr[j];
  74.             y[ihi] := yprr
  75.          END ELSE BEGIN
  76.             FOR i := 1 TO mpts DO
  77.                IF (i <> ilo) THEN BEGIN
  78.                   FOR j := 1 TO ndim DO BEGIN
  79.                      pr[j] := 0.5*(p[i,j]+p[ilo,j]);
  80.                      p[i,j] := pr[j]
  81.                   END;
  82.                   y[i] := func(pr)
  83.                END
  84.          END
  85.       END ELSE BEGIN
  86.          FOR j := 1 TO ndim DO p[ihi,j] := pr[j];
  87.          y[ihi] := ypr
  88.       END
  89.    END;
  90. 99:   END;
  91.